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Abstract: We compute the leading color, next-to-leading order QCD corrections to the 
dominant partonic channels for the production of a 1^ boson in association with three 
jets at the Tevatron and the LHC. This is the first application of generalized unitarity for 
realistic one-loop calculations. The method performs well in this non-trivial test and offers 
great promise for the future. 
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1. Introduction 

There are many multi-particle processes, knowledge of which through next-to-leading order 
(NLO) in QCD would be very desirable [1]. This statement, often repeated in the context 
of the forthcoming experiments at the LHC is, in fact, true even at the Tevatron. For 
example, the production rates for W{Z) + n jets, with n < 4 are well measured [2, 3] at 
the Tevatron^ but next-to-leading order QCD computations for such processes exist only 
forn < 2 [4]. 

Of course, there are good reasons for that. The NLO QCD computations for processes 
with a large number of external particles are difficult, both analytically and numerically; 
a list of well-known problems can be found in [1] . The need to overcome these difficulties 
has made the computation of one-loop multi-leg scattering amplitudes the focus of much 
research. In recent years, three main suggestions for possible solutions have emerged as a 
result. 

First, it was argued, and demonstrated by explicit computations, that traditional 
methods, where one starts from Feynman diagrams and proceeds through a Passarino- 
Veltman [5] style reduction, can be optimized and made highly efficient [6-16]. Second, it 
was also shown that pure numerical approaches to NLO computations are feasible [17-23]. 

The third idea is the use of generalized unitarity where one starts from on-shell tree- 
level scattering amplitudes and recycles them into loops. The idea of generalized unitarity 
was proposed in Ref. [24] more than ten years ago. Important physical results obtained 
using this method [25] have demonstrated both its potential and limitations. The tech- 
niques of applying generalized unitarity were significantly developed in recent years thanks 

^Currently, for total cross-sections, errors range from ten percent for W + 1 jet to fifty percent for W + 4 
jets. The error on VK + 3 jets production cross-section is about twenty percent [2, 3]. 
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TeV 


W+, LHC 


W~, LHC 


a [pb], /i = 40 GeV 
a [pb], = 80 GeV 
a [pb], /i = 160 GeV 


74.0 ± 0.2 
45.5 ± 0.1 
29.5 ± 0.1 


783.1 ± 2.7 
515.1 ± 1.1 
353.5 ± 0.8 


481.6 ± 1.4 

316.7 ± 0.7 
217.5 ± 0.5 



Table 1: The leading order total cross section for the production of a boson in association with 
three jets including both two quark and four quark processes vs. factorization and normalization 
scale. The results are obtained using the program MCFM. Cuts for the jets are px > 15 GeV, 
\r]\ < 2 at the Tevatron (^i = 1.96 TeV) and pr > 50 GeV, |r?| < 3 at the LHC (^i = 14 TeV). 
The CTEQ6L1 parton distributions which have as{Mz) — 0.13 are used. The quoted errors are 
statistical only. 

to important advances in Refs. [26-30] . These developments led to the design of two gen- 
eralized unitarity algorithms [31, 32]. These methods are seminumerical in the sense that 
they depend on the complete analytic knowledge of the relevant scalar integrals [33] . 

The computational algorithm suggested in Ref. [32] is employed in this paper; we will 
refer to it as D-dimensional generalized unitarity. Note that this method was recently 
used to obtain results not currently attainable with other methods, see e.g. Refs. [34-36]. 
However, an apparent weakness of generalized unitarity is that no result for a physical 
process has been obtained within this framework.^ This should be contrasted with the 
traditional tensor reduction approaches which never lost contact with phenomenology and 
are being constantly refined to accommodate new challenges. 

This is not a good situation for generalized unitarity which has to live up to the claim 
of its advocates that it is a more powerful method. The only way to address this potential 
criticism is to demonstrate the applicability of generalized unitarity in actual calculations 
of direct phenomenological interest, preferably in processes which are beyond the reach of 
traditional methods. We have chosen the production of a W boson in association with 
three jets for this purpose. The reasons for our choice are as follows: 

• the calculation of NLO QCD corrections to this process is of direct relevance since 
it is measured at the Tevatron [2, 3]; it is not possible to use the leading order (LO) 
prediction for a serious comparison of theoretical and experimental results because 
the LO cross section varies by as much as a factor of two under reasonable changes 
in renormalization and factorization scales, see e.g. Table |l]; 

• measurements at the Tevatron have shown that for W + n jets with n = 1 and 2, the 
data [2, 3] is well described by NLO QCD [4]; it is interesting to verify this also for 
three and higher numbers of jets; 

• W + 3 jet production is of interest for the LHC, being one of the backgrounds to 
model-independent searches for new physics using the jets plus missing energy signal; 

^We distinguish between generalized unitarity and application of the algorithm of Ref. [28] to Feynman 
diagrams. The latter method was employed for the computation of NLO QCD corrections to a relatively 
simple physical process pp —> VVV in [37]. 
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W,^/s = 1.96 TeV, 


= 14 TeV, 


^= 14 TeV, 


4f [pb] 


29.63 ± 0.04 


356.99 ± 0.77 


216.35 ± 0.40 


[pb] 


32.96 ± 0.02 


377.08 ± 0.79 


229.89 ± 0.42 


<^ [pb] 


16.13 ± 0.03 


147.60 ± 0.38 


94.91 ± 0.19 


< [pb] 


16.12 ± 0.02 


153.36 ± 0.36 


97.64 ± 0.22 



Table 2: Full and leading color cross sections for the production of a, W boson and three jets for 
the two- (2q) and four-quark (4q) processes, at leading order. Cuts, parton distributions and scale 
choices as in Tab. 0. The renormalization and factorization scale is set equal to 80 GeV. 

• the calculation of NLO QCD corrections to VF-l-3 jet production is highly non-trivial: 
there are 1583 Feynman diagrams including a significant number of high-rank six- 
point functions. There is no doubt that the computation of NLO QCD corrections 
to this processes is a challenging task for traditional diagrammatic approaches. 

In our opinion these reasons make -|- 3 jet production an ideal testing ground for 
the unitarity method and, recently, we made the first step in that direction. In Ref. [36] 
the three current authors, together with Giele and Kunszt, applied the generalized D- 
dimensional unitarity method to compute all one-loop matrix elements needed for the 
NLO corrections to -|- 3 jet production including two-quark (qqgggW) and four-quark 
(qqQQW) partonic processes [36]. Leading color two-quark amplitudes were also computed 
in Ref. [38]. 

With the virtual corrections in hand, two more steps are required to arrive at physical 
predictions fovW + 3 jet production. First, the virtual corrections should be integrated 
over the relevant phase-space. Second, we need to consider processes with one additional 
parton in the final state. When this parton becomes soft or collinear to other partons in 
the process, the final state that consists of four partons contributes to the final state with 
three hard jets. 

In spite of all the technical improvements described above, the computation of the 
matrix elements for TV -|- 5 partons at one-loop and Pi^ -|- 6 partons at tree-level and inte- 
grating them over the phase-space are very challenging tasks. For this reason it is useful to 
look for approximations, which help to reduce the technical complexity of the problem and 
are justifiable from a physics viewpoint. An obvious possibility is to consider the large- A'^c 
approximation. 

To check how well this approximation works, we study leading order results for l^-|-3 jet 
production. A compilation of results using the program MCFM [4] is given in Tables |||. 
It follows from these Tables that, both at the Tevatron and the LHC, two-quark processes 
dominate over four-quark processes.'^ At both colliders two-quark processes provide about 
70% of the observed cross-section, with four-quark processes being responsible for the re- 
maining 30%. Also, the large-A^r; approximation turns out to be good to about 10% for 
both the Tevatron and the LHC. We therefore conclude that a useful first step towards 

^We show numbers for fixed renormalization and factorization scales, but the relative decomposition 
into channels is largely scale- independent. 
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Tr , vs = 1.96 TeV, pp 


Vv^, a/s = 14 TeV, 


TUT— /~ 1 A T r 

W , vs = 14 TeV, pp 


r 1 1 

a [pb] 


32.96 ± 0.02 


377.08 ± 0.79 


229.89 lb 0.42 


99 [%] 


2.60 


6.80 


11.16 


95 [%] 


49.76 


37.90 


34.51 


59 [%] 


2.35 


37.86 


34.42 


^f; [0/1 


1 n on 


(Mi 


y.oo 


95 [%] 


3.63 


7.99 


9.38 


99 [%] 


0.0 


0.0 


0.0 


99 [%] 


0.0 


0.0 


0.0 


99 [%] 


21.52 


0.73 


0.60 


9-9 [%] 


0.26 


0.73 


0.60 



Table 3: The cross section for the production of a and three jets at leading color for the two 
quark processes and the percentages contributed by various incoming channels. Cuts and parton 
distributions as in Table |l|. The renormalization and factorization scale is set equal to 80 GeV. 



computing NLO QCD corrections to + 3 jet production cross-section is the calculation 
of those corrections for partonic processes with only two quarks in the initial and/or fi- 
nal state, in the large- A'^c approximation. Because the contribution of gg channel is small 
both at the Tevatron and the LHC, we may further limit ourselves to study processes with 
at least a quark or an anti-quark in the initial state, namely the six incoming channels 
qq,qq,qg,qg,gq, and gq. As we explain below, by working in the large- A'^c approximation 
and by considering the two-quark channels only, we can simplify the calculation signifi- 
cantly. 

The remainder of the paper is organized as follows. In Section |2| the computation of 
real emission corrections is described; the integration of the T^-|-6 parton matrix elements 
squared over the available phase-space is discussed in Section |3|. In Section |^ the calculation 
of the virtual corrections performed in Ref. [36] is reviewed. In Section |5| numerical results 
are presented. We conclude in Section ^ 

2. Tree-level processes and subtraction terms 

In this Section we discuss the computation of the relevant tree-level scattering amplitudes. 
We need these amplitudes to calculate the production cross-section for -|- 3 jets at tree 
level as well as the real emission correction to that process from the W + 4: parton final 
state. In what follows, we present matrix elements that describe the production of a 
boson, but everything that we say can be adapted to the case of W~ production, after 
obvious modifications. 

The scattering amplitude for the process 0^u + d + ng + can be decomposed 
into color-ordered amplitudes according to the equation 

<^^(l^,2,,3„..,n,) = 5"-' Yl iT--'''--T'''^'''X' Mlu,2d;a{3),,..,a{n)g). (2.1) 

a-£Sn-2 
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In Eq. ( |2.lD g is the strong coupling constant and Sn-2 denotes the (n — 2)! permutations 
of the gluons. Note that neither the W boson nor the electroweak coupUngs and CKM 
matrix elements are displayed in Eq. ( |2.lD . We employ a normalization of the color SU{3) 
generators such that Tr (T"T'') = 6""^. 

To calculate the production cross-section, we need to square the matrix element in 
Eq. ( |2.1| ) and sum over the color and spin degrees of freedom of the quarks and gluons. In 
the large- A'^c limit, individual color-ordered amplitudes do not interfere and we obtain the 
scattering amplitude squared 

Y,\A'ri'i-u,2d,3,..,n)\' = {g^y-'x^ ^ | ^(1^, 2,; ct(3), .., a(n))|2, (2.2) 

coljhel hel,S,i— 2 

where = {N^ — 1)N^^^. Note that we decided to keep some terms in the factor X„ 
which are subleading in the large- A^c limit. 

To arrive at the production cross-section, we need to choose a partonic initial state, 
square the scattering amplitude and integrate it over the phase-space available for the final 
state particles. Each choice of the initial state ij leads to a particular number of identical 
particles in the final state that we denote by Nij. When the integration over the phase- 
space available for final state particles is performed, we have to divide the result by the 
symmetry factor Sij = Nij\. If we combine the symmetry properties of the phase-space 
with the fact that no interference terms are present in Eq. ( |2.2D , we can reduce the number 
of scattering amplitudes that need to be calculated. 

For example, if the initial state is ud, extreme simplifications occur. In this case the 
final state is -|- (n — 2) g, and the number of identical gluons is N^^ = n — 2, leading to 
a symmetry factor S'^j = N^^ ! = (n — 2) ! = Sn-2- We therefore write 



a, 



ud 



coLhcl heLS, 



= {a'^T ^ / d</>giue^|^n(l«,2d;3g3,4g4,..,.ng„)|^, (2.3) 

hel 

where in the last step we used the symmetry of the (n — 2)-gluon phase-space, d(^giue, to 
argue that all color-ordered amplitudes give identical contributions to the cross-section. 
Since Eq. (2^) is the consequence of the fact that gluons are identical particles, it holds 
true independently of cuts or other restrictions imposed on partons in the final state. 

Other partonic channels can be simplified in a similar manner although, typically, we 
gain less compared to the ud initial state. For example, if we consider the initial state 
composed of a quark or anti-quark and a gluon, there are (n — 3) identical gluons in the 
final state. Therefore, the symmetry factor is Sqg = {n—3) ! and the number of independent 
color-ordered amplitudes that need to be considered is 8(^^-2)1 Sqg = (n — 2). For 3-parton 
final states n = 5, so that there are three independent amplitudes while for 4-parton final 
states, n = 6 and the number of independent amplitudes is four. When these numbers are 
compared to the ^3 = 6 and 5*4 = 24 independent amplitudes that would be required if 
fixed ordering of the final state particles were discarded, we see that the improvement is 
substantial. 
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When the matrix element with W + 4 partons in the final state is integrated over the 
phase-space subject to the requirement that three jets are observed, divergent results are 
obtained. These divergences arise from phase-space regions where one of the four partons in 
the final state becomes soft or two partons become collinear to each other; eventually, they 
cancel against similar divergences in the virtual corrections. To achieve this cancellation 
in practice, divergences in real emission contributions need to be extracted. To make 
those divergences manifest, we use a subtraction method [39] as formulated by Catani and 
Seymour [40] who proposed simple subtraction terms, which they called dipoles. However, 
we need to make minor modifications to the Catani-Seymour formalism because we work 
with amplitudes where the ordering of identical particles in the final state is fixed. 

To find dipole subtraction terms consistent with fixed ordering of the identical particles, 
it is convenient to follow the derivation in Ref. [40]. To this end, we calculate the soft limit 
of the amplitude squared, partial fraction the eikonal factors to isolate individual collinear 
limits, and extend the eikonal factors beyond the soft limits, taking Ref. [40] as an example. 
The dipoles that do not have soft singularities can be found by examining collinear limits 
of the contributing amplitudes. Although identical steps are required to find conventional 
dipoles, the difference between dipole terms that employ ordered and full amplitudes can 
be traced back to different soft limits and in the related necessity to go beyond the soft 
limit in a different way. 

We now illustrate the construction of the subtraction terms by considering the ud 
initial state. Upon ordering the gluons in the final state, the cross-section is determined 
by the square of a single color-ordered amplitude summed over helicities. In the actual 
computation of the cross-section, this amplitude squared is multiplied by the infrared-safe 
measurement function Fj that depends on the momenta of n partons. We therefore define 

P(2; 3, 4...n; 1) = \A„{U, 2^, 3, .., n)\^Fj{l, 2, 3, ...n), (2.4) 

hel 

where labels 1 and 2 denote incoming particles. Any gluon in the final state can become 
soft. We calculate the soft limit of the amplitude squared and obtain 

limP(2; 3, 4, 5, 6; 1) = s(5, 6, 1) V{2; 3, 4, 5; 1) + s(4, 5, 6) P(2; 3, 4, 6; 1) 

soft 

+s(3, 4, 5) V{2; 3, 5, 6; 1) + s{2, 3, 4)P(2; 4, 5, 6; 1). (2.5) 

The eikonal factor in Eq. ( |2.5D 

s{h3, k) = - — r (2.6) 

corresponds to the limit where momentum pj is soft.^ Performing partial fractioning and 
extending eikonal factors beyond the soft limit, we obtain the expression for the subtraction 
term 

Aub(2; 3, 4, 5, 6; 1) = 5 ® V{2; 3, 4, 5; 16) + Dll^<S) V{2; 3, 4, 56; 1) 



^We use a convention of treating all particles as if in the final state and as if all the momenta are 



outgoing. This allows us to write the soft limit in Eq. (2.5) in a symmetric way. 
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+i)f| 6 V{2- 3, 45, 6; 1) + -Dff 4 ® P(2; 3, 4, 56; 1) + ® V{2; 34, 5, 6; 1) 
+DII 3 ® P(2; 3, 4~5, 6; 1) + 5f| 2 » ^(2; 34, 5, 6; 1) + Z)f| 4 P(23; 4, 5, 6; 1). (2.7) 



The notation ij and j in Eq. (|2.7| ) are the standard notations, see Ref. [40] . The mapping 
between momenta p ^ p that is required to evaluate the right hand side in Eq. (2.7) is 
constructed in the same way as in Ref. [40]. The dipole functions Ln^' that we introduced 
in Eq. (2/7) are closely related to the original Catani-Seymour dipoles [40] and we explain 
the exact correspondence below. Before going into this, we point out that a modification 
of the subtraction terms is required to remove the symmetry between the emitter and 
emitted partons, inherent in final-final and final-initial dipoles in the original formulation 
by Catani and Seymour.^ In our notation, the non-integrable singular limit of the dipole 
Dij^k is associated with the soft limit of parton i whereas the soft limit of parton j does 
not introduce a non-integrable singularity. 

To make things clear, we give examples of dipoles that we employ in the present 
calculation. For final-final dipoles, where both emitter and emitted partons are gluons, we 
use 



Dfl^(g)V{2;..ij,k..-l) 



1 



KPiPj) 



1 



1 - 2:^(1 - yij^k) 



1 



2piPj 



A5i2;..ij^,~k.-l)Al{2;..ij„~k..;l), 



(2. 



where e = (4— D)/2 is the parameter of dimensional regularization, and yij^k = PiPj/{PiPj'^ 
PiPk +PjPk), Zj = PjPk/{PiPk +PjPk) and / = (1 - Zj)pi - zjpj. The momenta of particles 
ij and k are [40] 



Pi + Pj 



Vij.k 



1 - yij,k 



Pk, Pi 



Pk 



1 - yij,k 



(2.9) 



In Eq. (2^) A{..ij^...) denotes the tree level amplitude with the polarization vector for 
particle ij removed. 

For final-initial dipoles, where both emitter and emitted partons are gluons, we use 



Dfl^(2)V{..ij,..~a) 



1 



PiPjXija 



1 



1 Zj ~^ Xij^a) 



+(l-e] 



2piPj 



A^{..ij ....a)Al{...ij^, ...a), 



(2.10) 



where Xij^a = I + PiPj /{PiPa + PjPa) and Zj = PjPa/ {PiPa + PjPa) and / = (1 - Zj)pi - ZjPj. 
Note that the plus -sign between the first and the second term in the equation for Xij^a is 
the consequence of the all-outgoing momentum convention. 

For initial-final dipoles, where emitter and emitted parton are quark and gluon respec- 
tively, we use 



Dfl,^vCk..za)-. 



1 



'^{PiPa)Xia,k 



1 Xid h -\- U' 



(1 + e) - Xia,k{^ -e) 



|^5(^-ia)P,(2.11) 



^For initial-final and initial- initial dipoles this symmetry is not present to begin with and we use standard 
expressions for those dipoles. 
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where Xia^k = 1 + PiPk/ {PiPa + PkPa), Ui = PiPal {piPa + PkPa)- 

The subtraction terms for all other partonic channels are constructed along similar 
lines. Since more order ings contribute to the amplitude squared for partonic channels 
other than ud and du, more dipoles need to be considered to account for all the singular 
limits. Initial-initial dipoles appear for all channels except ud and du. Finally, let us note 
that the need to subtract certain dipoles cannot be established from the soft limit of the 
amplitude. In those case, the analysis of coUincar singularities of the amplitude squared is 
required to determine the dipoles that need to be subtracted. 

We also note that a simple but extremely useful modification of the Catani-Seymour 
dipoles was suggested by Nagy [41]. The idea is to limit the subtraction to a small region of 
phase-space available for final-state particles. To this end, final-final dipoles are multiplied 
hy 0{a — y), final- initial dipoles hy 9[a — [1 — x)), initial- final dipoles by 9{a ~ u) and 
initial-initial dipoles by 9{a — v) where y,x,u and v are standard variables used in [40]. 
This has the advantage that the subtraction is not performed if the kinematics of four- 
parton final state is far away from the singular limit, leading to a considerable saving in 
computing time, because the matrix element squared associated with the excluded dipole 
need not be computed. 

As we have seen, the construction of the dipoles relevant for our purposes is straight- 
forward and requires only small modifications compared to the original formalism of Catani 
and Seymour. The next step in the subtraction program - the integration of the subtracted 
terms over the unresolved phase-space - is even more straightforward since the integrals of 
our modified dipoles can be easily extracted from the results quoted in Ref. [40]. To see 
this, note that for initial-final and initial-initial dipoles, we do not introduce any modifi- 
cations relative to Ref. [40]. We do modify final-final and final-initial dipoles, but there is 
a simple relationship between our dipoles D and the ones in Ref. [40], -Dcs- For example, 
for final-final dipole, we can write 

Dcs{z,y)=D{z,y)+D{l-z,y). (2.12) 

To compute the integral of the subtraction term, we need to integrate D{z, y) over y and z 

1 

J dy dzf{z,y)D{z,y), (2.13) 


where f{z, y) is the weight function. The important property of the weight function is 
that it is symmetric with respect to z — > 1 — 2; transformation. Because of this symmetry 
property, we conclude that 

1 1 
2 y dy dz/(z, y)D{z, y) = J dy dzf{z, y)Dcs{z, y). (2.14) 
n 

Since the integral that appears on the right hand side of that equation is computed in [40], 
J D{z,y) can be easily extracted. A similar reasoning can be used to obtain integrals of 
the final-initial dipoles. 
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3. Integration over the phase-space 



The next issue to be discussed is the integration of the difference between the matrix ele- 
ment squared and the subtraction term over the entire phase-space allowed by the external 
cuts. Wc use VEGAS [42] to adapt the integration grid automatically but we still need to 
generate the parton kinematics carefully to ensure efficient sampling. 

In addition, when trying to integrate over the phase-space, we face a difficulty inherent 
in any subtraction method. To illustrate the issue, consider a matrix element squared that 
requires a subtraction of a particular final-final dipole to make it integrable. The dipole 
is described by two standard variables y and z. The matrix element squared has non- 
integrable singularities at y = 0, z = and at y = 0, z = 1; these singularities are 
removed by subtraction. The difference between the matrix element and the subtraction 
term is still singular, but these singularities are integrable. We assume that the difference 
scales as If-^/y and as l/s/z or — z. Although these are integrable singularities, 

in order to have the standard estimate of the integration error when using Monte-Carlo 
integration techniques, it is mandatory to change integration variables to absorb the square- 
root singularities into the measure. 

Since we have a left-over square-root singularity for every dipole that we need to 
subtract from the matrix element squared, we need to do a large number of variable trans- 
formations if we want to absorb all the singularities. Note also that different dipoles need 
to be subtracted for different initial states. This means that the required changes of vari- 
ables can not be done globally and, instead, we have to adopt a multichannel integration 
technique. To this end, for each partonic channel, we 

• randomly pick a dipole that contributes to a chosen channel - all dipoles are given 
equal weights; 

• generate the phase-space for the chosen dipole in such a way that the square-root 
singularity can be absorbed into the integration measure. 

More specifically, suppose we are interested in the contribution of a particular partonic 
channel to the production cross-section for which A^^ dipoles are required to make it finite. 
We need to compute 



where dcj)^ denotes the phase-space with four partons in the final state, \A\'^ is the matrix 
element squared and l^lg^btr subtraction term. We rewrite the integral as 




(3.1) 




(3.2) 



n=l 



where an are some constants and 




(3.3) 



m 
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In Eq. ( |3.3| ) J„ are Jacobian factors. In the actual calculation, we choose the coefficients 
an to be equal although, in principle, it is possible to optimize their choice iteratively. 

We associate each with a particular dipole contribution that needs to be subtracted 
from the matrix element squared. Suppose that n is a final-final dipole. We start by 
generating the phase-space for three massless partons and the W boson, assuming that the 
kinematics of the massless parton is characterized by a dpj_/p± distribution (for transverse 
momentum pT above the jet pT cut) and by a uniform distribution in rapidity.^ Then, we 
use the fact that the four-parton phase-space factorizes into the product of the three-parton 
phase-space and the dipole phase-space that is completely specified by three additional 
variables. We denote the momenta of the three-parton final state as p, and momenta of 
the four-parton final state as p. Additional variables needed to describe the kinematics of 
the 3^4 splitting are y, z and (p. We therefore write 

d(/>4(p) = d(/)3{pij,pk)dpi{pij,pk), (3.4) 

and 

dpi{pij,Pk) = ^ dzi dyij^k (1 - yij,k)- (3.5) 

Having generated the momenta of three-parton final state p, we use them to construct 
the momenta p according to the following formulae 

Pi = ZiPij + yij^kZjPk + k±, 

Pj = ^jPij + Vij^kZiPk - k±, (3.6) 
Pfc = (1 - yij,k)Pk, 

Pm^i,j,k ~ Pm- 



The transverse momentum reads kj_ = \k±\ (cos (p + sin cj) ) and | = y/yz{l — z)2pijpk- 
The two auxiliary vectors vi^2 are such that 2 = —1; V1V2 = 0, vi^2Pij = 0, fi,2Pfc = 0. 

The choice of the Jacobian factor J„ allows us to absorb the square-root singularities; 
for final/final dipoles we use the form suggested in Ref. [43]: 

J-i = (3.7) 

^yijjjj^^yziZj 

Note that this Jacobian can be written in terms of the four parton momenta p, since, 
as follows from Eq. (|3.7| ), we can express Zi, zj and yij^k in terms of these momenta, 

PiPk 1 PiPj /o 

- "--l-zi, yij,k = ■ — • (3.8) 



PiPk + PjPk PiPk + PjPk + PiPj 



This remark is important since, as follows from Eq. (3^), each integral In requires the 
knowledge of all Jacobians Jm including the ones with m ^ n. This, however, is not a 
problem because all those Jacobians are uniquely expressed through momenta p in the 
spirit of Eqs. (|3) and (|3l8|) . 



^This is the standard procedure in MCFM. 
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The Jacobian J„ can be absorbed into the measure by a simple change of variables 

<ipi{pij,Pk) PijPk X fon\ 
Jn " 2^ ^'^'^ ^ ~ ^ ' 

where 

^=^' y^3,k = ^4j^k, 2i = ^ (i - (1 - 2^i) Vi + 4^«(i - ^i)) , (3.10) 

with < fj,ij^k < 1 and < < 1. So, given random numbers for fj,,^, and 4> and the 
momenta of three-parton final state p, we can calculate z and y, the phase-space weight and 
the corresponding momenta p of the four-parton final state. We then calculate the matrix 
element squared as discussed in the previous Section and derive a particular contribution 
to the differential cross-section for + 3 jet production. 

Another case which requires comment involves dipoles where initial particles are present. 
We will consider the final-initial dipoles as an example. The general strategy is very similar 
to what has already been discussed. We start by generating momenta for the three-parton 
final state and the H^-boson and three random numbers z, x and (p. We denote the mo- 
menta of one of the partons in the final state by pij and the momentum of the spectator 
in the initial state by Pa- We then generate the additional parton momentum according to 
the following formulae 

. 1 — X _ 

Pi = ZPij + (1 - Z) Pa + k±, 

X 

Pj = (1 - z)pij + Z—^Pa - k±, 

Pa = -, (3.11) 

X 

Pm^i,a,k ~ Prn-i 

(3.12) 

where k^^^ = |A;_|_| (cos(/) -|- sin^ u^), = —1, vi^2Pij = 0, vi^2Pa = and = 

A/2^~p^j^(r^^^)(r^^. 

For the final-initial dipole, we employ the Jacobian 

Jn' = . (3.13) 

To absorb this Jacobian into the integration measure, we make a change of variables along 
the lines discussed in connection with final-final dipoles. Furthermore, we have to make sure 
that the new momenta of the initial parton pa does not exceed the momentum of the proton 
since, a priori, the variable x can assume any value between and 1. If the momentum of 
the initial state parton exceeds the momentum of the proton, the corresponding event is 
rejected. We deal with initial-final and initial-initial dipoles along similar lines. 

4. Virtual corrections 

A detailed description of the calculation of all one-loop amplitudes needed for the NLO 
correction to W^+3jcts at hadron colliders is given in Ref. [36]. Here we recall the elements 
of the discussion needed for the purpose of this paper. 
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At one-loop, using the color basis of Ref. [44] and neglecting contributions from closed 
fermion loops, the color decomposition for the process — > u + d+n g + W~^ can be written 
in terms of left primitive amplitudes [45] as 



A, 



l-loop , 



!«; 2(^, 3g, . . . , Tig) — g 



p=2 (TeSn-2 



JXlX2 



X (-l)M^(k, a{p)^, a{3)^, 2,, a{n)^, ...,a{p + 1) J 



.(4.1) 



In Eq. (^) for p = 2 the factor (T • • • T)^^*! becomes (T^^j^^i).^*! ^nd for p = n the factor 
{F ■ ■ ■ F)xix2 becomes 6x1x2- before neither the W boson nor the electroweak couplings 
and CKM matrix elements are displayed in Eq. ( [4.1D . In the leading color approximation 
we retain only the p = 2 term and Eq. (|4.1|) simplifies to 



yll-loOp/l 9, Q 



g, . . . ,ng 



9 



2:2 



x(-l)M^(l^,2rf,cT(n)^,...,a(3)^) 
The matrix element squared is then given by 

2 ^ ^ l-^n ^(fn; 2g(, 3g, . . . , ng,).A^ ' (Ifi, 2^^, 3g, .., 77,^)1 
coljhel 

= 2(iV2_i)^n-2(^2)n-i ^ 2,, 3„ . . . , n,)^;(l^, 2,, 3„ .. 

hel,Sn-2 



(4.2) 



(4.3) 



,ng)\ , 



where again we choose to keep some subleading color terms. We can now take advantage 
of the symmetry of the phase-space and fix the ordering of identical particles in the final 
state. The procedure is described in detail in Section p|. 



5. Results 

In this Section we present the results of our calculation of NLO QCD corrections to PF-l-3 jet 
production. We begin by describing computational aspects of the problem. The total 
time needed for the calculation of virtual corrections is determined by how many one-loop 
primitive amplitudes must be evaluated and how computationally expensive they are. At 
leading color, we only need to calculate the fastest primitive amplitude, for which about 
50 ms are required for a given momentum/helicity configuration.^ 

To compute the hadronic cross-section, we sum over six partonic channels. As ex- 
plained in previous Sections, we fix the ordering of gluons in the final state. Then, for 
partonic channels with quark and anti-quark in the initial state, we evaluate the virtual 
primitive amplitudes 2^ = 8 times, since we sum over two helicities of each of the three 
gluons. For the qg,gq,qg,gq channels we consider three different orderings of the final- 
state fermion, relative to the gluons. As the result, we evaluate 24 = 3 x 2'^ primitive 

^AU numbers are given for a computer with 2.33 GHz Pentium Xeon processor and Intel fortran compiler. 
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amplitudes per partonic channel with a single fermion in the initial state. Therefore, for 
the computation described in this paper, we need to compute the leading color primitive 
amplitude 112 times for each phase-space point. This translates into a total of 5.6 seconds 
per phase-space point. 

This time is too large to allow us to compute the virtual corrections on a dynami- 
cal, self-adapting grid. We therefore adopt the following strategy for computing virtual 
corrections. First, we compute the tree level cross-section with a large number of points, 
2 X 10^, and establish an integration grid.^ Once the grid is fixed, we compute virtual 
corrections by running three different jobs each with 10^ evaluations using different seeds 
to start VEGAS off. We then average the results of these three evaluations. With this 
procedure, Monte Carlo errors for virtual corrections are around 0.7 — 1%. 

For real corrections, we need about 10 msec per phase-space point to compute ma- 
trix element squared and the subtraction terms^. Since computation of the real emission 
correction is inexpensive, we calculate these corrections following the standard MCFM pro- 
cedure of first doing a pre-conditioning run and then a final run. We used 10 times 4 • 10^ 
points for the pre-conditioning run and five times 8 • 10^ points for the final run. With this 
number of events, the Monte Carlo integration errors for the real (subtracted) contribution 
are around 0.4 — 0.7%. 

As far as final results are concerned, errors coming from virtual and real corrections 
are comparable and the total Monte Carlo integration error is in the range of 1 — 4%; 
this is better than the theoretical uncertainty of the results estimated with a standard 
renormalization and factorization scale variation. 

We are now in position to describe the numerical results of the computation. However, 
before we enter into this discussion, we remind the reader that our results are approximate 
for the following reasons: 

• wc employ the large-A^c approximation to compute the scattering amplitudes; using 
the leading order cross-section as a guide, we estimate that this approximation is 
accurate to about 10 percent; 

• we include only the two-quark processes qqgggW and ignore the four-quark processes 
qqQQgW. Even within the two-quark processes, we do not consider the partonic 
channel with two gluons in the initial state. For the leading order cross-section, we 
find that the four-quark processes increase the cross-section by thirty percent so that 
omitting them gives results accurate to about thirty percent. 

Because of these approximations, we warn the reader that absolute results for cross- 
sections and distributions that we report below should be used with caution. We believe, 

®There are other ways to establish the grid. For instance, we can omit the computationally expensive 
parts of the virtual amplitudes and keep only logarithms of kinematic invariants that come from 1/e poles. 
We have checked that changing the strategy for establishing the grid has no bearing on the final result. 

^This number depends on the value of the parameter a that determines how often subtraction terms 
need to be calculated [41]. The quoted value corresponds to q = 0.01. 
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Figure 1: Inclusive W^^+S jet cross-section at the LHC and the i^-factor defined as = ctnlo /""lo 
as a function of the renormahzation and factorization scales. Jets are defined with kx algorithm 
with R = 0.7 and pr > 50 GeV. Jet rapidities satisfy Ir]] < 3. The LO and NLO cross-sections are 
computed with CTEQ6L1 and CTEQ6M parton distributions, respectively. 



however, that ratios of NLO and LO results for various observables are less sensitive to 
these omissions. 

The numerical results for VF-l-3 jet production at NLO are obtained using the CTEQ6m 
parton distributions [46] which have a value of as{Mz) = 0.118. The evolution of the 
coupling constant is performed using the two-loop beta function 

^/ N ,9/, ,/ N , 33 — 2n/- 153 — 19n/- , ^, 

«a,) = -ta|(l + f'"s). "--f^. '■'= 2»(33-2,v) - <'''' 

where, in the spirit of the large- A'^c approximation, we set the number of light flavors n/ 
equal to zero. The kx jet algorithm with R = \J /S.cfP' + Ar/^ = 0.7 and px > 15 GeV {pT > 
50 GeV) is used to define jet cross sections at the Tevatron and the LHC, respectively. We 
employ default MCFM choice for electroweak parameters and the CKM matrix elements; 
they can be found in Ref. [4]. 

In Figs. |l[P we present total cross-sections and JC- factors, defined as K = unlo/o'lO) 
for 1^ -|- 3 jet production at the LHC and the Tevatron as a function of the factorization 
and the renormalization scales which we set equal to each other //^ = fip = fx. At the 
LHC, the NLO cross-section shows remarkable independence of the scale /i, unlike the LO 
result. The equality of LO and NLO cross-sections occurs at fiQ ~ 160 GeV. Because the 
dependence of the LO cross-section on the unphysical scale is strong, the NLO corrections 
are typically large. For example, choosing ^ = mw to compute the LO cross-section for 
-|- 3 jet production at the LHC, leads to NLO QCD corrections of the order of —50%. 



- 14 - 



Tevatron 




40 80 120 160 200 240 
1^ 



Figure 2: The inclusive + 3 jet cross-section at the Tevatron and the ii'-factor defined as 
K = (Tnlo /o'lo a-s a function of the renormalization and factorization scales /x. Jets are defined 
with kx algorithm with R = 0.7 and pr > 15 GeV. Jet rapidities satisfy \ri\ < 2. The LO and NLO 
cross-sections are computed with CTEQ6L1 and CTEQ6M parton distributions, respectively. 



For the Tevatron, the situation is different. First, the dependence of the NLO cross- 
section on the renormalization and factorization scales is sizeable although it is significantly 
reduced compared to the leading order cross-section. In addition, as follows from Fig. |2| the 
equality of leading and next-to-leading order cross-sections occurs at a scale ^ 50 GeV 
which is much smaller than the LHC case discussed above. This is not unexpected since 
both the center of mass energy and the p± cut for jets is smaller at the Tevatron which 
leads to a much softer spectrum of jets compared to the LHC case. 

It is interesting to note that gross features of NLO QCD corrections to -|- 3 jet pro- 
duction, such as scales at which leading and next-to-leading order cross-sections coincide, 
are very similar to what was observed in NLO QCD computation of -|- 2 jets [4, 47]. 
What differs between two and three jet production is the price one pays for making an 
infelicitous choice of scale in the LO result. Because the dependence on /x of crj^^g is 
stronger than of cry^_^2 jet' -f'^-factors for 1^ -|- 3 jet decrease or increase stronger when one 
moves away from ^ = ^q. 

Finally, we present selected results for differential distributions at the LHC. We choose 
the renormalization and factorization scales to be 160 GeV since this minimizes the inclusive 
X-factor. In Fig. |3| we plot the distribution in the variable Ht defined as the sum of 
transverse energies of jets, the missing transverse energy and the transverse energy of the 

lepton Ht = ^ E^j + ^^j"^"" + El. The variable Ht measures the overall hardness of a 

j 

particular event and can be employed in model-independent searches for New Physics. As 
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Figure 3: The distribution of the transverse energy Ht = J2 ^^-j + E™^'^^ + and the K factor 

3 

defined a.s K = (dcr^^°/diJT) /{da^'^/dHT) in W+ + 3 jet inclusive production at the LHC at 
leading and next-to-leading order. Rcnormalization and factorization scales are set to 160 GeV. 



illustrated Fig. ^ the i/T-distribution becomes softer at NLO, at least in comparison to 
the leading order result calculated at fixed scale fi = 160 GeV. 

In Fig. 1^ we present the transverse momentum distribution of the third hardest jet in 
the inclusive production of -|- 3 jets at the LHC. In the range of p± shown in the plot, 
the shapes of p± distributions at leading and next-to-leading order are nearly identical. 

6. Conclusions 

In this paper, we apply the method of generalized D-dimensional unitarity [32] to compute 
NLO QCD corrections to -|- 3 jet production at the Tevatron and the LHC. There 
are two reasons that make this result an important benchmark in the field of one-loop 
QCD computations for hadron collider physics. First, this is the only application of the 
idea of generalized unitarity in a fully realistic one-loop computation that goes beyond 
calculation of one-loop helicity amplitudes at a fixed point in phase-space. Second, our 
result is one of the very few computations of one-loop corrections to six-parton processes 
at hadron colliders - the current research frontier in NLO QCD. It is remarkable that 
the method achieved this benchmark without a problem; this assures us that generalized 
unitarity is a practical computation method that can be applied to other, perhaps even 
more complicated, processes. 

Looking into the near future, we expect to refine our computation in two ways. First, 
we expect to include the four-quark partonic channels in the large- A'^c approximation; this 
is an important step for realistic phenomenology. Further down the road, we may want to 
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Figure 4: The transverse momentum distribution of tiie third hardest jet in the inclusive + 3 
jet production at the LHC. The renormalization and factorization scales are set to 160 GeV. 



extend the computation beyond the leading color approximation. Estimating the increase 
in computer running time required to go beyond the large Nc limit, we find that about 
two minutes per point will be needed to compute virtual corrections to the matrix element 
squared. At face value, this is feasible, but computationally expensive. However, one can 
imagine various improvements, including Monte Carlo sampling over helicities and colors, 
that should lead to an appreciable improvement in the speed of the program. 

Finally, we would like to say a few words about phenomenology. Since we did not 
consider the four-quark channels in this paper, we decided not to pursue very detailed phe- 
nomenological studies. However, the numerical results that we do report are instructive 
since they give an idea about potential significance of NLO QCD effects in + 3 jet pro- 
duction at the Tevatron and the LHC. Our computation shows that NLO QCD effects are 
large and can reach ±50%, if unfortunate, but not unreasonable, choices of the renormal- 
ization and factorizations scales are made in a computation based on leading order matrix 
elements. Note that the probability that an unfortunate scale choice is made increases for 
a larger number of jets since the production cross-section at LO becomes a steeper function 
of the renormalization and factorization scales. The only way to cure this problem is by 
computing NLO QCD corrections. For processes like pp ^ W + 4: jets or pp ^ tt + 2j 
this will be complicated no matter what method is used, but we believe that generalized 
unitarity will be up to the task. 
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